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Propagation of radio waves in the earth-ionosphere waveguide is considered for the case where the 
lower boundary is an inhomogeneous smooth surface. An integral equation for the problem is formu- 
lated in a direct fashion by utilizing the compensation theorem. After some simplifications, several 
special cases are considered explicitly. For example, in the case of a two-section path consisting of 
a long stretch of sea and a short section of land, a relatively simple working formula is obtained. The 
result shows that the modal excitation factors at VLF for an all sea path are significantly reduced when 
the foreground is poorly conducting. Another special case considered is when the propagation path 
is all sea except for a short intermediate land section. In this case, it is found that energy from low- 
order modes will be transferred to high-order modes with a subsequent reduction of field strength. 



1. Introduction 

A great deal of attention has been given to the prob- 
lem of predicting groundwave fields for mixed land/sea 
paths. Furthermore, extensive calculations have been 
made which show the interdependence of the various 
parameters. Since a recent review of progress in 
this area is now available [Wait, 1964], it is not neces- 
sary to discuss this particular topic here. However, 
it is rather surprising that, in propagation via iono- 
spheric reflections, little attention has been given to the 
influence of an inhomogeneous earth. It is the pur- 
pose of this paper to consider this problem with special 
reference to VLF radio propagation. 



2. Formulation 

The mutual impedance z m between two vertical 
antennas at A and B (separated by a great circle dis- 
tance d) located on a spherical earth of radius a is 
considered. In order to account for the presence of 
the ionosphere, an equivalent reflecting layer is lo- 
cated at height h. For present purposes, we assume 
this layer is characterized by a surface impedance Z\ 
which does not vary along the path. If over the earth's 
surface between the terminals A and B, the surface 
impedance is Z everywhere, the mutual impedance 
may be expressed as a sum of modes as follows 



Zm — Zq 



d)a 



sin (did) 



1/2 



W{Z, d), 
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where z is the mutual impedance between the dipoles 
A and B if they were located on a perfectly conducting 
flat ground plane and separated by a distance d. In 



the above expression, W(Z, d) is an attenuation func- 
tion defined by [Wait, 1962]: 



W(Z,d) = 



(h/k) l 



nlA V exp 



(ikdf 



An, (2) 



where C n is the cosine of a complex angle, and A w is 
the excitation factor for modes of order n. The CVs 
are solutions of a modal equation which involves the 
following dimensionless parameters: 

{-t n ) ll2 = C n {kal2yl* 

y = kh(kal2)-^ 

q=-i(ZI>no)(kal2yl a 

qi =-i{Zilr)o){kal2yi\ 



where k = 27r/wavelength, and 170— 12077 ohms. 
After making a number of simplifying assumptions, 
the modal equation mentioned above may be written 
in the form [Wait, 1962] 



W' 2 {tn) — qW 2 (tn) 
w[(t n ) — qWi(t n ) 



W'l(tn — Jo) + qjWijtn — Jo) 
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where w\{t) and w 2 (t) are Airy functions while the 
primes indicate derivatives with respect to the argu- 
ments. An alternative form of (2) is 



ro n 



ixtn) A„ 



(4) 
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where 



X = (dla)(kal2yi*. 



Numerical values of t n and the excitation factor 
A n are now available for a wide range of parameters 
[Wait and Spies, 1964]. Thus, for the purposes of 
this communication, the solution of the waveguide 
problem with constant wall impedances Z* and Z is 
taken to be known. What is of interest here is the 
extension to the same waveguide when the surface 
impedance of the lower boundary is a function of 
position along the path connecting A and B. For 
example, as indicated in figure 1, it is now assumed 
that, over some area S on the earth's surface, the sur- 
face impedance is Z r which may be different from Z. 
The method of approach is very similar to that used by 
the author [Wait, 1964] in studying mixed-path ground- 
wave propagation. As in that case, the problem may 
be formulated in terms of the compensation theorem 
in the form derived by Monteath [1951]. For example, 
the mutual impedance z' m between dipoles A and B 
over an inhomogeneous ground of variable surface 
impedance Z' is given by 



m 



z 'm — z m+n I \(Z' —Z)H at -H' M dS 



(5) 



where z m is the mutual impedance if the surface of the 
earth were homogeneous with surface impedance Z 
everywhere. The tangential magnetic field of dipole 
A over the homogeneous earth is H a t while the tan- 
gential magnetic field over the inhomogeneous earth is 
H^. The currents in the dipoles are both taken equal 
to /o for convenience. The surface of integration S 
extends over the region of the earth which is charac- 
terized by a surface where Z' differs from Z as indi- 
cated in figure 1. 

The formal equivalence of (5) with the formula of 
mixed-path groundwave theory is a result of the im- 
pedance boundary condition on the upper wall of the 
waveguide. In other words, the surface impedance 
Z\ is assumed to be the same for both the homogeneous 
and for the inhomogeneous earth cases. If this were 
not permissible, there would be an additional surface 
integral over the upper boundary of the waveguide. 




Figure 1. 



Plan view of the inhomogeneous region on the earth's 
surface. 



To simplify the present problem, the mutual impedance 
z' m is defined in terms of an attenuation W'{Z, Z', d) 
such that 



z m~ z 
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sin (d/ 
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W'{Z, Z', d) 



(6) 



in analogy to (1) for the homogeneous earth. The 
next step is to express the tangential magnetic field 
vectors, at the point of integration P in terms of at- 
tenuation functions. Thus 



and 



ikloha e _ ik8 / 1+ 1 



2tts 



ikloh 
2wl 



b e - ik i 



1 + 



iks 



iki) 



s/a 



sin (s/a) 



l/a 



1/2 



F(Z, S) (i n Xi s ) 

(7) 



1/2 



sin (l/a) 

JT(Z,Z\Z)(i n xi,), (8) 



where s and / are great circle distances from A and B 
to P and i s and i* are unit vectors in the directions of 
increasing s and Z, respectively. In the above, h a 
and hb are the effective heights of the dipoles A and B 
while i n is a unit vector normal to and into the surface 
of the spherical earth. 

Equation (5), when combined with (6), (7), and (8), 
leads to the equation 



W'(Z,Z' 9 d) = W(Z 9 d)+ e 



ikd 



-I I Gis, 

o J J s 



g-ik(s+l-d) 



^ttTt] J J s ~ si 

X (Z f - Z) W(Z, s)W'(Z, Z', cos 8 dS, (9) 
where 8 is the angle subtended by i s and i* and where 
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This is a two dimensional integral equation for the 
unknown function W(Z, Z', d). To solve such an 
equation directly appears to be hopeless. Therefore, 
some simplifications are made at this stage in order 
to achieve tractability. 



3. Approximate Form of the Integral 
Equation 

We note that the function exp [— ik(s + 1 — d] is 
rapidly varying compared with other factors in the 
integrand. Therefore, one may expect that the prin- 
cipal contribution to the integrand will occur when 
s + l — d, provided that the surface impedance con- 
trast Z'—Z does not change rapidly in a direction 
transverse to the path. Therefore, in the other factors 
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in the integrand, / may be replaced by a, and s may 
be replaced by d — a, where a is the great circle dis- 
tance from B to the point Q on the great circle between 
A and B. (The arc QP is perpendicular to AB.) 
Furthermore, over most of the range of integration, 
(1/ks) and (1/kl) may be neglected compared with 
unity and similarly cos 8 may be replaced by — 1. 
The latter approximations are obviously violated when 
the terminals are near discontinuities of the surface 
impedance Z\ 

The reduced form of the integral equation is 

ikd 
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where 
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where all quantities, except the exponential factor, 
vary with a only. The exponent s + l — d is now ex- 
pressed in terms of the angle coordinates fl and (3 
with reference to figure 2; these are defined by 

a = al a = QB/a and /3 = QP/a 
From spherical trigonometry 

cos (BOP) = cos 11 cos 13 
and 

cos (AOP) = cos (6-fy cos (3. 
Thus, 



BOP=-=H+ 
a 



cot n 



/3 2 + terms in /3 4 , (3 6 



(14) 



and 



AOP=-=#-n+ 

a 



cot(O-il) 



+ terms in /3 4 , /3 6 , 
Therefore, to a first order in /3 2 , 



(15) 
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[cot n+ cot (0-n)]/3 2 . 
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Figure 2. Coordinates for describing the integration over the 
spherical surface of the earth. 



For the present development it is further assumed 
that Z' — Z does not vary in the transverse direction 
(i.e., with (3) over the surface S. Thus, Z' is regarded 
only as a function of a. 

Within the region of validity of the present approxi- 
mations, the element of area dS may be approximated 
by a 2 d£ld(3. The integral now has the form 



r'(Z, Z\d) = W(Z, d) 

f« 2 Z \a)-Z 
(d-a) 



ikda f™ Z_ 
27D7o J«, a 

p 2 («) 
X 

J/3,(a) 
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exp [-i(kal2) [cot (1 + cot {6-n)](3 2 ]df3da, 

(17) 

where /3i(a) = y\(a)la, /32(a) = y2(a)la. In this case, 
the limits of the surface S are y\(a) ^ y ^ y<i(a) and 
a\ =^ a ^ «2 as indicated in figure 1. After a change 
of the /3 variable, (17) may be written in the form 

W\Z,Z\d) = W(Z, d) - f^V /2 d \ a2 Z ' {a) ~ Z f(a) 

„ r(Z,rf-q)tT(Z,Z',tt) m 

r(uu ui)da 



a(d — a) 



I -) + cot 
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where 



F(ai, ud^dlW* f" 2 exp (-i7ru 2 /2)du, 
ui = (£a/7r) 1/2 [cot H + cot (0 - n)] 1/2 /3i, 



(18) 

(19) 
(20) 



and 



u 2 = (ka/7T)^ 2 [cot a + cot (6 - ft)] 1/2 /3 2 . (21) 



As indicated before, the rapidly varying function 

exp [— ik(s + l — d)] in the integrand of (9), deter- 

(16) mines the portions of the earth's surface which are 

significant. The phenomenon may be described in 
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terms of Fresnel zones. These are determined by the 
locus of the points where 

k{s + l-d) = mTTl2{ior m = l, 2, 3 . . .) 

or, to within a good approximation, 

(Aa/2)[cot (a/a) + cot ((d-a)la)]p 2 = mTrl2. 

The width of the first Fresnel zone at any point a is 
then obtained from 

fm = 2pa = (2maky' 2 [cot (a/a) + cot ((d-cO/a)]" 1 ' 2 . 

A 

The maximum width, denoted f m , occurs where a==d/2. 
Explicitly, 



where 



f m = (mkdl2yiWI 2 , 



y _ tan [rf/(2a)] [rf/(2a)]« 

A ~ <f/(2a) _1+ 6 + - * * * 



(22) 



(23) 



Under most conditions Z may be replaced by unity 
even when d is comparable with the earth's radius a. 
Within this approximation, the Fresnel zones are el- 
lipses and the semiminor axes are/ m /2, while the semi- 
major axes are (a 7 /2) + (wX/8). In the classical sense, 
the "first Fresnel zone" corresponds to m = 2. It is 
immediately evident that F(u\, ^2) being replaced by 
unity in (18), is equivalent to saying that the width of 
the surface S extends to several Fresnel zones on both 
sides of the propagation path. 

It is convenient to rewrite (18) in the following form: 

w\z, z\ d) = w{z, d)-(^) 112 



x 



i a2 Z'(a)- 
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F(u u u 2 ) 
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(24) 



where we have made use of the trigonometric identity 



sin (d/a) 
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(25) 



When F(ui, 112) is set equal to one, (24) bears a formal 
equivalence to the one dimensional integral equation 
developed for groundwave propagation over mixed 
paths [e.g., Wait, 1964]. It is interesting to note 
that (24), in the form given, is not restricted to dis- 
tance d such that d/a«l. This comes about be- 
cause of the normalization factor [(d/a)/ sin (did)] 112 
which is not included in the definition of the attenua- 



tion function W(Z, d). In the case of groundwave 
propagation, this distinction is not of any consequence 
since d/a is small in any case. However, in VLF 
propagation in the earth-ionospheric waveguide, the 
ratio d/a may be comparable with unity and the nor- 
malization factor mentioned above may exceed unity 
by a significant amount. 

Equation (24) is in a reasonably tractable form for 
direct numerical calculation of the attenuation func- 
tion W\Z, Z', d) when the limits of the surface S and 
the surface impedance function Z' are specified. For 
present purposes, some simple limiting cases will be 
considered rather than attempting a frontal assault 
on (24). 

4. Two Section Path 

An important special case of the general mixed path 
problem is when Z f is sectionally homogeneous, as 
indicated in figure 3. For example, if the path be- 
tween A and B is characterized by a constant surface 
impedance Z, from to d— d\ and Z\ from d— d\ to d, 
the integral equation (24) simplifies to 



W\Z, Z u d) = W(Z 9 d}- 



ikd\ 1,2 (Zi 
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when d\ > and when the function F(ui, 112) has been 
replaced by unity. Of course, if d\ < such that A and 
B are both over the surface of the earth of surface 
impedance Z, it is seen from (24) that 

W'{Z,Z U d) = W(Z, d). 

At least this is true to within the stationary phase 
approximation which, in effect, reduces the area inte- 
gration to a line integration. As a result of this reason- 
ing it is equally justified to replace W'(Z, Zi, a) where 
it occurs in the integrand of (26) with W(Z\, a) which 
is the attenuation function for propagation from the 
point B to distance a over a homogeneous earth of 
surface impedance Z\. 

Equation (26) with the simplification indicated in the 
preceding paragraph may be written in terms of the 
dimensionless coordinates mentioned earlier. Thus, 



EQUIVALENT REFLECTING LAYER 




Figure 3. Sectional view of earth-ionosphere waveguide for a two- 
section path, such that Zi # Z. 
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for the two section path, 

W'{q, q uX ) = W{q, x) + Q0"* (<?, ~ q) 

y [*' W{q,X-xW(.qu-k) 
Jo [(X~X)X] 

where 



1/2 



dx, (27) 



q = -i{Zly }0 ){kal2YI\ Ql =- J(Z,/T, )(A:a/2) 1/3 , 
X = {dla)(kaliyi\ Xi = (dila)(ka/2yl 3 , and 
X = {ala)(kal2yi\ 
Because of (4), we may write 



and 



' (g,X-X) 2 /77V/ 2 v 

(x-x)" 2 -y«W £,,/" 
(x) 1/2 yo\i/ ^,, 2 , 3 .. 



e -'(x-x« K 



-W„ 



(28) 



(29) 



where A n is the excitation factor corresponding to the 
earth of surface impedance Z while A m is the excita- 
tion factor corresponding to the earth of surface 
impedance Z\. The t n are roots of (3) while T m are 
roots of an equation identical to (3) if q is replaced 
by q x . 

On inserting (28) and (29) into (27), the integration 
with respect to x ma Y De readily carried out to yield 

W{q,q u x) = W(q, X ) 

+ 3 e««l*(7r X y' 2 (q, ~ «) 2 2 AJLr** ., r , • 

(30) 

The physical significance of this result is best seen by 
examining the special case where the distance d— d\ 
or x~ Xi i s sufficiently large that only the n=\ term 
is needed. Thus 



W\q,q l ,x)^y o (j) m Aie- i x<i 



+ ± e i3*l4 {7rx y,2 {qi - q)Aie - txtl 

Ai — " =~ h A2 " 



The first term on the right of (31) is the attenuation 
function for a path which is homogeneous throughout 
its length (with a surface impedance Z). The remain- 
ing terms, proportional to q\ — q or Z\—Z are correc- 
tions which result from the inhomogeneity extending 
over the path of length d\ (proportional to x0- The 
terms proportional to A 2 , A 3 . . . etc., represent 
conversion of energy in the waveguide from mode 1 
to mode 2, 3 . . . etc. 

A somewhat simpler approach to the two section 
problem is appropriate when the section of the path 
of length d\ is very small compared with the total length 
of d (i.e., Xi <<: X)- Then, from (27), it is seen that 

W'(q,quX)^W(q,x) + J^3i 



X 



\l l W{q, X -x)W{quX)^ (32) 



which is a slight simplification. When the latter re- 
sult is applied to a long sea path (of length d — d\) and 
a short land path (of length d\), further simplifications 
are possible. For example, W{q, \~ X)> m the inte- 
grand, is replaced by W{q, x) ana \ thus, 



r'< 9 ,g,,x) = ^fa,x)[i + «] 



where 



a= 



2L 



Utt) 1 ' 2 



/; 



W(q u x) 

(X)" 2 



dx- 



(33) 



(34) 



On inserting the mode series expansion given by (29) 
into (34), the integration with respect to x is readily 
carried out. However, because of the assumed small- 
ness of Xu the resulting expansion would be very poorly 
convergent. An alternative approach is to recognize 
that W(x, Q\) ma Y be replaced by the groundwave at- 
tenuation function for x ranging from to Xi- 

The form appropriate for short distances is [Wait, 
1964] 



JT(quX)= 2 A m e im ^qf(x) ml \ (35) 



m=0, 1,2... 



where 



A = 1, A x =- i V^rT A t =-2, A 3 = i V^ (l +^A 



3l 29? 



—J, • • • etc. 



(ti-ti) 



(ti-t 2 ) 






Using (35), the integration indicated in (34) is carried 
(**■) out to yield the following series form for the correction 
factor: 
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where 



and 



p}/ 2 =(*d 1 /2) 1 / 2 (Z 1 /7 ?0 )e-^ 4 



(36) 



(37) 



qi =-i(kal2yi*(Ziho)- 

It is seen from this series expansion that the earth 
curvature only influences the correction factor ft 
through the higher order terms in powers of p\ /2 . 
Also, if Z and Z\ correspond to sea and land respec- 
tively, the factor 1 — (Z/Zi) may be replaced by unity. 
Furthermore, at VLF, displacement currents in the 
ground are negligible, which means that (j gx + ie gi 0) 
may be replaced by a gi in the definition of p\. Thus 

pP~(llk)(dil*gWnm*l* 

is a real quantity. Then, to a first order 



^here 



ft' --pi and n r, = -(2l7T 1 ' 2 )p 1 J 2 . 



Thus, the fractional reduction in the amplitude of the 
field is p\ while the phase lag is increased by 2(pi/7r) 1/2 
radians. For example, if the length d\ of the land 
path is 100 km, for a land conductivity cr 9l of 1 milli- 
mho/m and a frequency of 15 kc/s (i.e., A — 20 km), it 
is easily found that 

ft' --2.62 X 10- 2 and 0"=- 0.129 rad = - 7.4°. 

For this example, the short land section has a negligi- 
ble effect on the propagation over the total distance d. 

5. Propagation Across a Strip 

An interesting situation occurs when the path 
between A and B is homogeneous with surface im- 
pedance Z everywhere except for a relatively short 
stretch of length d 2 where the surface impedance is 
a constant Z 2 as indicated in figure 4. Assuming 
that the inhomogeneity is effectively a strip of in- 
finite transverse dimension, F(ii\, u 2 ) in (24) may be 
replaced by unity. The resulting integral equation 
for the attenuation function may then be written 



W' = W(q, X ) + 






(<72-<?) 



-W{q,X-xW'{q,q2,X) 



[(x-x)x? /2 



d X , (38) 




Figure 4. Sectional view of the earth-ionosphere waveguide for an 
intermediate section where Z 2 7^ Z. 

where q 2 =-i(ka/2)^ (Z 2 /r ?0 ), Xi = (ka/2)^ d x \a and 
X2 — (ka/2yt 3 d 2 /a. The function W\ as it occurs in 
the integrand in the preceding equation, may be 
identified as the attenuation function appropriate for 
propagation from the point B to a variable point on 
the strip (i.e., \i +X2 > X > Xi)- I n accordance with 
the previous discussion, this particular function W 
may be replaced by the appropriate form for a two 
section path which, in effect, ignores reflection at 
the boundary x = X2 (i.e., at distance di + d 2 from B). 
If the series representation given by (30) is used, the 
subsequent integration leads to a triply infinite series 
for the resultant attenuation function for the path 
A to B. Convergence of this type of expansion is 
satisfactory provided that the distance parameters 
X, Xu and X2 are all somewhat greater than one. 

A somewhat simpler approach to the strip problem 
is to regard the midsection as a perturbation to 
the homogeneous path. In this approximation, 
W(Q,„Q2, x) in the integrand of (38) is replaced by 
Jf(<7» X\ which is the attenuation function for propaga- 
tion over a homogeneous earth of surface impedance 
Z. Furthermore, if the strip is relatively narrow 
(i.e., X2 < Xi and x)» the integrand may be replaced 
by its value at the midpoint of the strip. Thus, 



W'^W(q, X )+ A 

\1TI 



where 



1/2 
(Q2 - q) 



X2 



[Xo(*-Xo)] 1/2 



*lT(q,X-XoW(q 9 Xo), (39) 



Xo = Xi + (X2/2) = (ka/2Y*(do/a), do = dt + (A/2). 

The perturbation term, which is proportional to q 2 — q, 
may now be regarded as a first order or single scatter- 
ing from the strip. 

Using the modal representation of the form given by 
(28), it readily follows from (39) that 



where 



Ar-^(i} /2 (,2-,)X2 



(40) 



xX- 



-i(x ~ XOH 



"An 2 e-^mA m . (41) 
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If x is sufficiently large, only the n = \ term is needed 
and the above representation for AW may be written 
in the more meaningful form, 

-^r-- — (<? 2 -g)x2[Ai + A2e- J,( '2-<i)xo 

+ A 3 e~ i( '3 ~ ' l)x ° + Aie-m - h)xo + . . . ] (42) 



or, what is the same thing, 

AW_ (Z 2 - 
W - \ V o 



where 
and 



^\y[A 1 + A 2 e-^-^i)^o 

+ A 3 c-^3-5i)*do+. . .], (43) 



s m = a-ci l yi*^i-(ci l /2), 
c m =(-t m y l2 {2ikayi\ 



Thus, as indicated, the strip will modify the strength of 
the first mode by a factor proportional to Ai and, at the 
same time, it will produce higher order modes propor- 
tional to A2, A 3 , . . . etc. 

For negligible displacements in the earth, the multi- 
plicative factor in (43) may be written 

Z<l — Z\ di 



Vo 



n-m'^-miP"'- <«> 



where cr 92 is the conductivity of the earth over the 
strip of width d 2 . For €r g >a g2 , ^2=100 km, /i = 70 
km, /= 0)l2u = 15 kc/s, cr g2 = \ milli-mhos/m, it follows 
that 



AW 
W 



^-0.040 e i7r[4 






(45) 



For this example, the excitation factors Ai, A 2 , . . . are 
not appreciably different from unity and they are 
nearly real. Thus the strip produces approximately a 
0.03 fractional diminution of the amplitude of first 
mode and a change of phase of the order 0.03 rad. 
The higher modes are excited with a relative strength 
of about 0.04. The resultant effect of the higher 
modes, of course, depends on the magnitude of the 
electrical distance kdo and the relative phase velocities 



of the modes. In general, the higher modes are attenu- 
ated as a result of the increasing value of the imaginary 
part of S m as m increases. As indicated above, 
numerical values of the excitation factors A /M and the 
propagation factors S m are available [Wait and Spies, 
1964] for a variety of conditions appropriate in the 
VLF range. 

6. Concluding Remarks 

It would appear that the inhomogeneity of the ground 
is an important factor in the propagation of VLF radio 
waves in the earth-ionosphere waveguide. In certain 
practical applications, such as to navigation systems 
and worldwide communications, the influence of 
inhomogeneous land sections on the path may alter 
significantly the normal behavior of the transmission. 
In particular, abrupt changes in conductivity from 
sea to land may convert appreciable amounts of energy 
to the higher modes. Not only will this change the 
resultant attenuation, but the effective phase velocity 
will be modified as a result of modal interference. 



I thank K. P. Spies and R. L. Gallawa for their 
helpful comments. 
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